function [EZ_welfare_PHI,V1,p_grid] = EZ_welfare_PHI_v4(Wi,CAPH,MU,maxreps,tolFP,rho,psi,alpha,invG,G,prob_start)T = size(Wi,1);Nstates = size(MU,1);Pspot_N = PGR_prices_v2(MU,CAPH,rho,T,tolFP,maxreps);% prices can only fall% so max price comes from the first period% i will assume no one starts in the states that would yield negative% consumption so i avoid issues with the CES in the EZ preferences maxP= max(Pspot_N(:,1).*(Pspot_N(:,1)<Wi(1)));step = round(maxP/5000);step = max(step,1);p_grid = [0:step:maxP];%maxP = max(max(MU));%step = round(maxP/1000);%p_grid = [0:step:maxP];P_grid = repmat(p_grid,Nstates,1);%mortalityP_grid(Nstates,:)=0;ngrid = size(P_grid,2);%PHI premiums    P_bar_f= Pspot_N(:,T);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% period T-1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t  = T-1;  % C_f = C{t+1}|c_t, lambda_{t+1}% premium grid(cols) and future lambda in rowsP_f = P_grid.*(P_grid<=repmat(P_bar_f,1,ngrid))+repmat(P_bar_f,1,ngrid).*(P_grid>repmat(P_bar_f,1,ngrid));Y_f = repmat(Wi(T),Nstates,1);% mortalityY_f(Nstates,1) = 0;Y = repmat(Wi(t),Nstates,1);Y(Nstates,1) = 0;C_f = Y_f - P_f;% The value function in each future lambda (lambda_t+1)% For terminal period is just (1-rho)*C_f% this is F(c_t,0)V_f = (1-rho)^(1/(1-1/psi))*C_f;%V_f = CES_agg(C_f,0,psi,rho);% Take expectations over lambda capt = CAPH(:,:,t);GV_f = G(V_f,alpha);GV_f(Nstates,:)=0;EGV_f = capt*GV_f;  %inverse operatorinvGEGV_f = invG(EGV_f,alpha); invGEGV_f(Nstates,:)=0;% for our given p; and for each state; what is V% using CES aggregatorV = CES_agg(Y-P_grid,invGEGV_f,psi,rho);Vfirst = V;t= T-2;while (t>=1)  %t  %if (mean(mean(imag(V)))>0)  %    break  %end  capt = CAPH(:,:,t);    P_bar_f = Pspot_N(:,t+1) ;  Y = repmat(Wi(t),Nstates,1);  Y(Nstates,1) = 0;    %% for every possible premium/lambda in the future;   %% this is V in the future  V_f = V;  % If you land in the state in the row (j); and   % you currently have premium (i)  % your future price is P_f(i,j)  % so this is P_{t+1} | lambda_{t+1}, P_t  P_f = P_grid.*(P_grid<=repmat(P_bar_f,1,ngrid))+...    repmat(P_bar_f,1,ngrid).*(P_grid>repmat(P_bar_f,1,ngrid));  % round to the next integer  % this makes the whole thing an approximation    P_f = round(P_f/step)*step;  Col_P_f = P_f/step+1;    %now I need to read V(P_{t+1}}|lambda_{t+1},c_t}    % if future state is in row r  % and future price is equal to P  % then I need to read Vf at:   % row r and   % the column that corresponds to P in the P_grid   % which means I need to read the following index        index =(Col_P_f-1)*Nstates+repmat([1:Nstates]',1,ngrid);    GV_f = G(V_f(index),alpha);  GV_f(Nstates,:)=0;  EGV_f = capt*GV_f;    invGEGV_f = invG(EGV_f,alpha);   invGEGV_f(Nstates,:)=0;  V = CES_agg(Y-P_grid,invGEGV_f,psi,rho);   t = t-1;  end  %starting price  P1 = Pspot_N(:,1);  % this is V(State,pgrid)  delta = (repmat(p_grid,Nstates,1) -repmat(P1,1,ngrid)).^2;  mindelta = min(delta,[],2);    P1_loc = delta==repmat(mindelta,1,ngrid);  V1 = sum(V.*P1_loc,2);  % CERTAINTY EQUIVALENT  % TAKING INTO ACCOUNT THE FIRST-PERIOD UNCERTAINTY    CAPC =eye(length(MU(:,:,1)));   psurv = zeros(T,1);  psurv(1,1) = 1;  d = zeros(T,1);  d(1,1)=1;  for t = 2:T      CAPC = CAPC*CAPH(:,:,t-1);      psurv(t,1) = 1-CAPC(1,length(MU(:,:,1)));      d(t,1) = psurv(t,1)*rho^(t-1);  end    GV1 = G(V1,alpha);  EGV1 = prob_start*GV1;  V0 = invG(EGV1,alpha);  a = 1/(1-rho)*V0^(1-1/psi)/sum(d);  EZ_welfare_PHI= a^(1/(1-1/psi));    %  a = G(V0,alpha)/(1-rho)/sum(d);  %  EZ_welfare_PHI_v3= invG(a,alpha)end